function cfdInitialShockTube
%--------------------------------------------------------------------------
%  Written by LiuHaHa @ NWPU, 2022 01 
%  Contact me at: liuzhikan@mail.nwpu.edu.cn
%==========================================================================
% Routine Description:
%   This function initialize the flow field using shock tube 
%--------------------------------------------------------------------------
global cfdFlow slvpara;

theNumberOfElements = cfdFlow.mesh.numberOfElements;
elementCentroids = cfdFlow.mesh.elementCentroids;

if(slvpara.ND == 2)
    for i=1:theNumberOfElements
        if(elementCentroids(i,1) < slvpara.initPara(1))
            cfdFlow.Phy(i,1) = slvpara.initPara(2);
            cfdFlow.Phy(i,2) = slvpara.initPara(3);
            cfdFlow.Phy(i,3) = slvpara.initPara(4);
            pressure = slvpara.initPara(5);
            cfdFlow.Phy(i,4) = pressure/(cfdFlow.gama-1) + 0.5*(cfdFlow.Phy(i,2)*cfdFlow.Phy(i,2) + ...
                cfdFlow.Phy(i,3)*cfdFlow.Phy(i,3))*cfdFlow.Phy(i,1);
        else
            cfdFlow.Phy(i,1) = slvpara.initPara(6);
            cfdFlow.Phy(i,2) = slvpara.initPara(7);
            cfdFlow.Phy(i,3) = slvpara.initPara(8);
            pressure = slvpara.initPara(9);
            cfdFlow.Phy(i,4) = pressure/(cfdFlow.gama-1) + 0.5*(cfdFlow.Phy(i,2)*cfdFlow.Phy(i,2) + ...
                cfdFlow.Phy(i,3)*cfdFlow.Phy(i,3) )*cfdFlow.Phy(i,1);
        end
    end
else
    for i=1:theNumberOfElements
        if(elementCentroids(i,1) - slvpara.initPara(1) < 1e-15)
            cfdFlow.Phy(i,1) = slvpara.initPara(2);
            cfdFlow.Phy(i,2) = slvpara.initPara(3);
            cfdFlow.Phy(i,3) = slvpara.initPara(4);
            cfdFlow.Phy(i,4) = slvpara.initPara(5);
            pressure = slvpara.initPara(6);
            cfdFlow.Phy(i,5) = pressure/(cfdFlow.gama-1) + 0.5*(cfdFlow.Phy(i,2)*cfdFlow.Phy(i,2) + ...
                cfdFlow.Phy(i,3)*cfdFlow.Phy(i,3) + cfdFlow.Phy(i,4)*cfdFlow.Phy(i,4))*cfdFlow.Phy(i,1);
        else
            cfdFlow.Phy(i,1) = slvpara.initPara(7);
            cfdFlow.Phy(i,2) = slvpara.initPara(8);
            cfdFlow.Phy(i,3) = slvpara.initPara(9);
            cfdFlow.Phy(i,4) = slvpara.initPara(10);
            pressure = slvpara.initPara(11);
            cfdFlow.Phy(i,5) = pressure/(cfdFlow.gama-1) + 0.5*(cfdFlow.Phy(i,2)*cfdFlow.Phy(i,2) + ...
                cfdFlow.Phy(i,3)*cfdFlow.Phy(i,3) + cfdFlow.Phy(i,4)*cfdFlow.Phy(i,4))*cfdFlow.Phy(i,1);
        end
    end
end





























